Analisi dell'innalzamento del livello dei mari

Davide Giardini - Progetto di "Big Data in Geographic Information Systems", anno accademico 2022-2023

L'argomento dell'innalzamento del livello dei mari è da tempo affrontato dall'opinione pubblica, parzialmente a causa del sempre maggiore interessamento verso le tematiche legate al riscaldamento globale, ma in particolare anche per i potenziali effetti che questo fenomeno potrebbe avere sulle aree costiere. Le cause che più spesso vengono attribuite all'innalzamento della superficie degli oceani sono lo scoglimento dei ghiacciai e l'espansione dell'acqua dovuta all'aumento della sua temperatura.
Andremo quindi a esplorare tre diversi modelli: uno sul livello della superficie del mare, uno sulla temperatura dell'acqua e un'ultimo sull'area coperta dai ghiacciai. Lo scopo finale sarà quello di valutare la correlazione e, possibilmente, la causazione delle ultime due variabili sulla prima.

Operazioni preliminari sui modelli singoli

# Import libraries
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns

import cartopy.crs as ccrs

import xarray as xr
import nc_time_axis

import numpy as np
import pandas as pd

Modello 1 - Sea Conservative Temperature

La temperatura conservativa è una proprietà termodinamica dell'acqua di mare. È derivata dall'entalpia potenziale ed è raccomandata secondo lo standard TEOS-10 in sostituzione della temperatura potenziale poiché rappresenta in modo più accurato il contenuto di calore nell'oceano.

Per questa variabile utilizziamo il modello "CMIP6.CMIP.CNRM-CERFACS.CNRM-CM6-1.historical" disponibile qui.
Il modello è stato costruito a seguito del progetto: "Voldoire, Aurore (2018). CMIP6 simulations of the CNRM-CERFACS based on CNRM-CM6-1 model for CMIP experiment historical".

# Apro ed esploro il modello 1
modfile1 = 'C:/Users/david/Desktop/GeoData/Project/bigthetaoga_Omon_CNRM-CM6-1_historical_r17i1p1f2_gn_185001-201412.nc'
d1 = xr.open_dataset(modfile1)
d1

<xarray.Dataset>
Dimensions:      (time: 1980, axis_nbounds: 2)
Coordinates:
    sector       |S255 b'global'
  * time         (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:...
Dimensions without coordinates: axis_nbounds
Data variables:
    time_bounds  (time, axis_nbounds) datetime64[ns] 1850-01-01 ... 2015-01-01
    bigthetaoga  (time) float32 3.954 3.957 3.958 3.958 ... 4.039 4.041 4.044
Attributes: (12/51)
    Conventions:            CF-1.7 CMIP-6.2
    creation_date:          2019-07-15T09:30:04Z
    description:            CMIP6 historical
    title:                  CNRM-CM6-1 model output prepared for CMIP6 / CMIP...
    activity_id:            CMIP
    contact:                [email protected]
    ...                     ...
    xios_commit:            1442-shuffle
    nemo_gelato_commit:     49095b3accd5d4c_6524fe19b00467a
    arpege_minor_version:   6.3.2
    parent_variant_label:   r1i1p1f2
    history:                none
    tracking_id:            hdl:21.14100/84a180f0-104e-4149-aa29-1c2bc876b254

Analisi della Serie Storica mensile

Per poter decomporre la serie storica mensile nelle sue tre componenti (trend, stagionalità e residui) dobbiamo apportare delle modifiche alla struttura dei dati in modo che possano essere passati alla funzione "seasonal_decompose" della libreria "statsmodels".

# Trasformiamo i dati in dataframe
ts_bigthetaoga = d1.drop_vars('time_bounds').to_dataframe()['bigthetaoga'].to_frame()

# Aggiungiamo un indice temporale identico alla variabile "time" già presente nel modello, ma che possa essere interpretato
# dalla funzione di decomposizione
ts_bigthetaoga['month'] = pd.to_datetime(pd.date_range('1850-01', '2015-01', freq='M'))
ts_bigthetaoga = ts_bigthetaoga.set_index('month')

La funzione seasonal_decompose ci permette di decidere tra i due tipi di decomposizione (additivo e multiplicativo) tramite l'opzione "model". La scelta del modello è determinata dall'analisi dei residui: cerchiamo il modello che risulti in residui con range minore, varianza costante, e che non presentino pattern.

from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse

# Additive Decomposition 
bigthetaoga_decom = seasonal_decompose(ts_bigthetaoga.bigthetaoga, model='additive', extrapolate_trend='freq')

# Plot
plt.rcParams.update({'figure.figsize': (10, 10)})
bigthetaoga_decom.plot()
plt.show()

Ci concentriamo ora sulla stagionalità annuale:

plt.rcParams.update({'figure.figsize': (9, 5)})
bigthetaoga_decom.seasonal.groupby(bigthetaoga_decom.seasonal.index.month).mean().plot()
plt.title("Seasonal Component of Global Average Sea Water Conservative Temperature")
plt.show()

La temperatura conservativa media globale del mare presenta una ben determinata componente stagionale. Questa presenta un picco in corrispondenza del mese di Aprile e un minimo in corrispondenza di quello di Settembre.
La variazione però è relativamente contenuta: 0.012 gradi attorno al trend.

Rimozione della Componente Stagionale

Andiamo ora a isolare il trend, rimuovendo quindi la componente stagionale, calcolando una media annuale della temperatura.

d1= d1.resample(time="Y").mean()

plt.rcParams.update({'figure.figsize': (9, 5)})
d1.bigthetaoga.plot()
plt.title("Global Average Sea Water Conservative Temperature")
plt.show()

Il grafico risultante ci dimostra un trend decisamente crescente, che è cresciuto, dal 1850 al 2014 di 0.09 gradi. Il grafico sembra anche mostrare un'accellerazione di questo trend dagli anni 1960.

Modello 2 - Sea Surface Height Above Geoid

Per valutare l'innalzamento del livello dei mari utilizziamo la variabile "zos": "Sea Surface Height Above Geoid". Come dice il nome, la variabile calcola la distanza verticale della superficie oceaninca da una superficie di riferimento, che in questo caso è il geoide. Il geoide è una superficie di geopotenziale costante, con la quale il livello del mare dovrebbe coincirdere se gli oceani fossero a "riposo". Fonte: NVS.

Per questa variabile utilizziamo il modello "zos_Omon_BCC-CSM2-MR_esm-hist" disponibile qui.
Il modello è stato costruito a seguito del progetto: "Wu, Tongwen; Chu, Min; et al. (2018). BCC BCC-CSM2MR model output prepared for CMIP6 CMIP esm-hist.".

# Apro ed esploro il modello 2
modfile2 = 'C:/Users/david/Desktop/GeoData/Project/zos_Omon_BCC-CSM2-MR_esm-hist_r1i1p1f1_gn_185001-201412.nc'
d2 = xr.open_dataset(modfile2)
d2

<xarray.Dataset>
Dimensions:    (time: 1980, bnds: 2, lat: 232, lon: 360)
Coordinates:
  * time       (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
  * lat        (lat) float64 -81.5 -80.5 -79.5 -78.5 ... 86.5 87.5 88.5 89.5
  * lon        (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    latitude   (lat, lon) float32 -81.5 -81.5 -81.5 -81.5 ... 88.26 88.05 87.83
    longitude  (lat, lon) float32 0.5 1.5 2.5 3.5 ... 60.95 63.36 65.25 66.77
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (time, bnds) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00
    lat_bnds   (lat, bnds) float64 -81.5 -81.0 -81.0 -80.0 ... 89.0 89.0 89.5
    lon_bnds   (lon, bnds) float64 0.0 1.0 1.0 2.0 ... 358.0 359.0 359.0 360.0
    zos        (time, lat, lon) float32 ...
Attributes: (12/49)
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    branch_method:          Standard
    branch_time_in_child:   0.0
    branch_time_in_parent:  2145.0
    comment:                The model integration starts from the esm-piContr...
    ...                     ...
    title:                  BCC-CSM2-MR output prepared for CMIP6
    tracking_id:            hdl:21.14100/d97c7956-1712-409f-a7ca-58d4b42c929b
    variable_id:            zos
    variant_label:          r1i1p1f1
    license:                CMIP6 model data produced by BCC is licensed unde...
    cmor_version:           3.3.2

Analisi Geografica

Ci concentriamo prima su un'analisi geografica. Andiamo a plottare l'altezza media del mare sopra al geoide.

# Situazione media
plt.rcParams.update({'figure.figsize': (8, 5)})
p0 = d2.zos.mean(dim='time').plot(transform=ccrs.PlateCarree(),
                            subplot_kws={'projection': ccrs.Robinson()})
plt.title("Average Sea Surface Height Above Geoid through time")
p0.axes.coastlines()
plt.show()

Dal grafico risultante notiamo come l'area in cui l'altezza media del mare tende ad essere superiore a quella del geoide si addensa intorno all'equatore, mentre ai poli l'altezza media del mare tende ad essere inferiore.

Andiamo ora a verificare se sono avvenuti cambiamenti durante gli anni: costruiamo una mappa che mostri le differenze tra i primi 20 anni di osservazioni (1950-1979) agli ultimi anni (1995-2005). Le aree in rosso andranno a indicare i punti in cui il livello del mare è aumentato maggiormente, al contrario quelle in blu corrisponderanno con i punti in cui il livello del mare è diminuito.

# Analisi delle differenze tra i primi 20 anni e gli ultimi 20 anni
plt.rcParams.update({'figure.figsize': (8, 5)})
Anom = d2.sel(time=slice('1995-01-01', '2015-01-01')).mean(dim='time') - d2.sel(time=slice('1850-01-01', '1870-01-01')).mean(dim='time')
p0 = Anom.zos.plot(transform=ccrs.PlateCarree(), subplot_kws={'projection': ccrs.Robinson()})
plt.title("Anomalies of Sea Surface Height Above Geoid \n between first and last 20 years")
p0.axes.coastlines()
plt.show()

Il grafico risultante è quasi per la totalità di colore rosso, dimostrando un aumento generale del livello dei mari.
L'unica area che ha dimostrato un abbassamento dei livelli del mare è quella nelle prossimità dell'Antartide. Al contrario, a latitudini appena superiori a questa fascia di colore blu, risiede la zona corrispondente al maggiore innalzamento dei livelli dei mari (constrassegnata dal rosso più intenso).

Dall'analisi grafica notiamo anche che il modello sembra presentare degli artefatti, soprattutto in relazione del Polo Nord. Le coordinate descritte nel modello non sembrano combaciare bene con quelle delle funzioni di plot. Vari tentantativi per risolvere il problema sono stati fatti, ma senza successo.
Considerando il nostro interesse verte su una serie storica della media globale della variabile stessa, e non su una meticolosa analisi geografica, abbiamo deciso di utilizzare comunque il modello.
Il problema si ripeterà con il modello successivo.

# Usiamo uno "slice" delle latitudini per permetterci di concentrarsi sul polo Nord
p1a = d2.zos.isel(time=0,lat=slice(210,233)).plot(transform=ccrs.PlateCarree(),
                subplot_kws={'projection': ccrs.Orthographic(central_longitude=0.0, central_latitude=90)},
                                                            cmap='Blues')
p1a.axes.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x1aa3e0032b0>

Costruzione Serie Storica Bilanciata

Per costruire la serie storica dobbiamo calcolare la media globale dell'altezza del livello del mare sopra al geoide. Per fare ciò dobbiamo ovviamente calcolare la media attraverso ogni cella. Dobbiamo però prima definire dei pesi in modo da rendere conto della diminuzione dell'area delle celle all'avvicinarsi verso i poli. In questo modo potremo calcolare una media pesata e raggiungere una serie storica più veritiera.

fig = plt.figure(figsize=(9,5))  
d2.zos.mean(axis=1).mean(axis=1).plot()
plt.title("Sea Surface Height Above Geoid prima del bilanciamento")
plt.ylabel("Sea Surface Height Above Geoid [m]")
plt.show()

# Costruzione di variabili utile rappresentanti le dimensioni dei dati
nmonths = len(d2.time)
nlat = len(d2.coords['lat'].values)
nlon = len(d2.coords['lon'].values)

# Definizione dei pesi in modo da tener conto della diminuzione dell'area coperta dalle celle all'avvicinarsi dei poli
weights = np.tile(np.cos(np.deg2rad(nlat)), (nlat,360))

# Costruzione di un array che possa contenere i dati mensili del modello
data_glob = np.empty([1, nmonths],dtype=float)
data_glob[:,:] = np.nan

# Costruzione di una maschera per rimuovere i valori nan (dati dalle celle corrispondenti a terre ferme)
masked_data = np.ma.masked_array(d2.isel(time=0).zos, np.isnan(d2.isel(time=0).zos))

# Loop sui mesi e calcolo della media pesata per il modello 2
for t in range(0,int(nmonths)):  
    masked_data = np.ma.masked_array(d2.isel(time=t).zos, np.isnan(d2.isel(time=0).zos))
    data_glob[0,t] = np.ma.average(masked_data,weights=weights)

# Costruzione della visualizzazione per il modello 2
fig = plt.figure(figsize=(9,5))  
plt.plot(pd.date_range('1850-01', '2015-01', freq='M'),data_glob[0])
plt.title("Sea Surface Height Above Geoid dopo del bilanciamento")
plt.ylabel("Sea Surface Height Above Geoid [m]")
plt.show()

Analisi della Serie Storica mensile

# Costruzione di un dataset idoneo all'utilizzo della funzione "seasonal_decompose"
ts_zos = pd.DataFrame(data_glob[0], columns=["zos"])
ts_zos['month'] = pd.to_datetime(pd.date_range('1850-01', '2015-01', freq='M'))
ts_zos = ts_zos.set_index('month')

# Additive Decomposition 
zos_decom = seasonal_decompose(ts_zos, model='additive', extrapolate_trend='freq')

# Plot
plt.rcParams.update({'figure.figsize': (10, 10)})
zos_decom.plot()
plt.show()

Ci concentriamo ancora una volta sulla componente stagionale mensile.

plt.rcParams.update({'figure.figsize': (9, 5)})
zos_decom.seasonal.groupby(zos_decom.seasonal.index.month).mean().plot()
plt.title("Seasonal Component of Global Average Sea Surface Height Above Geoid")
plt.ylabel("Sea Surface Height Above Geoid [m]")
plt.show()

La serie stagionale dell'altezza del livello dei mari presenta un andamento quasi completamente opposto di quella della temperatura conservativa. Essa infatti presenta un picco di 0.0015 nel mese di Settembre ed un minimo di -0.0022 ad Aprile.

Rimozione della Componente Stagionale

zos_year = [sum(data_glob[0][i:i+12])/12 for i in range(0,len(data_glob[0]),12)]

ts_zos = pd.DataFrame(zos_year, columns=["zos"])
ts_zos['year'] = pd.to_datetime(pd.date_range('1850-01', '2015-01', freq='Y'))
ts_zos = ts_zos.set_index('year')

plt.rcParams.update({'figure.figsize': (9, 5)})
ts_zos.plot(legend=False)
plt.title("Global Average Sea Surface Height Above Geoid")
plt.ylabel("Sea Surface Height Above Geoid [m]")
plt.show()

Dal precedente grafico, ottenuto applicando al modello una media annuale della variabile di riferimento, un trend positivo lineare che ha portato il livello del mare al di sopra del geoide da un minimo di 2.22 cm nel 1850 ad un massimo di 9.74 cm nel 2015.

Modello 3 - Ice Area

Per quanto riguarda lo scoglimento dei ghiacchiai, utilizziamo la variabile "simass", cioè la "massa totale di ghiaccio marino diviso per l'area della cella". Considerando che utilizzeremo, come proxy dello scoglimento del ghiaccio, la massa di ghiaccio totale, ci aspettiamo di avere correlazioni opposte a quanto sarebbe stato intuibile utiizzando una variabile che misurasse lo scoglimento. Ci aspettiamo insomma una correlazione negativa di questa variabile con le altre, considerando con una diminuzione dell'area del ghiaccio corrisponde a un maggiore scoglimento dello stesso.
Da notare anche che la variabile è già divisa per l'area della cella, per cui non consideriamo necessaria un'operazione di pesaggio nel momento del calcolo della media globale, come invece era stato per il modello 2.

Per questa variabile utilizziamo il modello "BCC-CSM2-MR_esm-hist" disponibile qui allo stello link del precedente.
Il modello è stato costruito a seguito dello stesso progetto del modello 2: "Wu, Tongwen; Chu, Min; et al. (2018). BCC BCC-CSM2MR model output prepared for CMIP6 CMIP esm-hist.".

# Apro il modello 3
modfile3 = 'C:/Users/david/Desktop/GeoData/Project/simass_SImon_BCC-CSM2-MR_esm-hist_r1i1p1f1_gn_185001-201412.nc'
d3 = xr.open_dataset(modfile3)
d3

<xarray.Dataset>
Dimensions:    (time: 1980, bnds: 2, lat: 232, lon: 360)
Coordinates:
  * time       (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
  * lat        (lat) float64 -81.5 -80.5 -79.5 -78.5 ... 86.5 87.5 88.5 89.5
  * lon        (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    latitude   (lat, lon) float32 -81.5 -81.5 -81.5 -81.5 ... 88.26 88.05 87.83
    longitude  (lat, lon) float32 0.5 1.5 2.5 3.5 ... 60.95 63.36 65.25 66.77
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (time, bnds) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00
    lat_bnds   (lat, bnds) float64 -81.5 -81.0 -81.0 -80.0 ... 89.0 89.0 89.5
    lon_bnds   (lon, bnds) float64 0.0 1.0 1.0 2.0 ... 358.0 359.0 359.0 360.0
    simass     (time, lat, lon) float32 ...
Attributes: (12/49)
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    branch_method:          Standard
    branch_time_in_child:   0.0
    branch_time_in_parent:  2145.0
    comment:                The model integration starts from the esm-piContr...
    ...                     ...
    title:                  BCC-CSM2-MR output prepared for CMIP6
    tracking_id:            hdl:21.14100/04f7ca22-6892-418e-8dae-e386dea2dbc3
    variable_id:            simass
    variant_label:          r1i1p1f1
    license:                CMIP6 model data produced by BCC is licensed unde...
    cmor_version:           3.3.2

Analisi Geografica

Per prima cosa effettuiamo un'analisi geografica della distribuzione della massa di ghiaccio. Senza alcuna sorpresa, notiamo che la maggior parte di essa è concentrata nel polo Nord, mentre un'altra parte si addensa intorno al contidente dell'Antartide.
Una cosa importante da fare per questa analisi geografica è porre un limite massimo allo span dei colori utilizzando l'opzione "vmax" della funzione "plot". Questo perché la distribuzione della massa di ghiaccio è fortemente assimmetrica: pochissime celle raggiungono valori molto alti. Se non usassimo questa opzione, lo span della legenda coprirebbe l'intero range dei valori, e così quasi tutto il grafico risulterebbe bianco.

plt.rcParams.update({'figure.figsize': (8, 5)})
p0 = d3.simass.mean(dim='time').plot(transform=ccrs.PlateCarree(), 
                                 subplot_kws={'projection': ccrs.Robinson()}, 
                                 cmap='Blues', vmax=1000)
plt.title("Average Total Mass of Sea Ice through time")
p0.axes.coastlines()
plt.show()

Come seconda cosa, invece, calcoliamo le anomalie tra i primi 20 anni e gli ultimi 20 anni di osservazioni, per poter analizzare i cambiamenti della distribuzione delle masse di ghiaccio. Il grafico sottostante, perciò, mostrerà in rosso le celle in cui l'area del ghiaccio è diminuita, e in blu le celle in cui essa è aumentata.

# Analisi delle differenze tra i primi 20 anni e gli ultimi 20 anni
plt.rcParams.update({'figure.figsize': (8, 5)})
Anom = d3.sel(time=slice('1995-01-01', '2015-01-01')).mean(dim='time') - d3.sel(time=slice('1850-01-01', '1870-01-01')).mean(dim='time')
p0 = Anom.simass.plot(transform=ccrs.PlateCarree(), subplot_kws={'projection': ccrs.Robinson()}, vmax=200, cmap = 'RdBu')
plt.title("Anomalies of Total Mass of Sea Ice \n between first and last 20 years")
p0.axes.coastlines()
plt.show()

Come ci aspettavamo, il grafico prende una colorazione quasi del tutto rossa, a indicazione di una generale diminuzione delle aree coperte dai ghiacciai globalmente. Ovviamente, le aree rosse si addensano ai due poli, in quanto, come abbiamo visto prima, sono le zone con più presenza di ghiacciai, e che hanno quindi subito una diminuzione maggiore.

Concentrandosi sull'area Antartica è interessante vedere come la differenza tra la quantità di ghiaccio presente nel mare sia notabile semplicemente plottando le due situazioni una a fianco all'altra, anche senza costruire un grafico che disegni le differenze tra i due periodi. I due grafici rappresentano la situazioni media dei primi e degli ultimi 10 anni.
È evidente una generale diminuzione dell'area dei ghiacciai nella regione attorno al continente antartico nell'ultimo periodo.

d3_sliced1 = d3.simass.isel(lat=slice(0,30)).sel(time=slice('1850-01-01', '1860-01-01')).mean(dim='time')
d3_sliced2 = d3.simass.isel(lat=slice(0,30)).sel(time=slice('2005-01-01', '2015-01-01')).mean(dim='time')

p1a = d3_sliced1.plot(transform=ccrs.PlateCarree(),
                      subplot_kws={'projection': ccrs.Orthographic(central_longitude=0.0, central_latitude=-90)},
                      cmap='Blues', vmax=500)
p1a.axes.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x1aa3ebe83d0>

p2a = d3_sliced2.plot(transform=ccrs.PlateCarree(),
                      subplot_kws={'projection': ccrs.Orthographic(central_longitude=0.0, central_latitude=-90)},
                      cmap='Blues', vmax=500)
p2a.axes.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x1aa3ea70580>

Analisi della Serie Storica mensile

ts_simass = d3.simass.mean(axis=1).mean(axis=1).to_dataframe()

ts_simass['month'] = pd.to_datetime(pd.date_range('1850-01', '2015-01', freq='M'))
ts_simass = ts_simass.set_index('month')

# Multiplicative Decomposition 
simass_decom = seasonal_decompose(ts_simass.simass, model='multiplicative', extrapolate_trend='freq')

# Plot
plt.rcParams.update({'figure.figsize': (10, 10)})
simass_decom.plot()
plt.show()

plt.rcParams.update({'figure.figsize': (9, 5)})
simass_decom.seasonal.groupby(simass_decom.seasonal.index.month).mean().plot()
plt.title("Seasonal Component of Sea-Ice mass per Area")
plt.show()

La serie stagionale della massa di ghiaccio presenta invece un andamento del tutto simile a quella della temperatura conservativa. Essa infatti presenta un picco di 1.364 kg/m2 nel mese di Maggio ed un minimo di 0.658 kg/m2 a Settembre.

Rimozione della componente Stagionale

d3= d3.resample(time="Y").mean()

plt.rcParams.update({'figure.figsize': (9, 5)})
d3.simass.mean(axis=1).mean(axis=1).plot()
plt.title("Global Average Sea Ice Area")
plt.show()

Fra le tre variabili presentate, la massa di ghiaccio è quella che presenta un trend meno netto se si considera l'intero span temporale. Quello che è molto visibile, però, è un trend decrescende partito dagli anni 1970 e ancora presente negli ultimi anni di osservazioni.

Confronto tra i 3 modelli

# Costruzione di un dataset che contenga tutte le variabili
ts_zos['bigthetaoga'] = d1.bigthetaoga.to_dataframe()['bigthetaoga']
ts_zos['simass'] = d3.simass.mean(axis=1).mean(axis=1)
ts_zos.head()

	zos	bigthetaoga	simass
year
1850-12-31	0.025079	3.952904	105.652702
1851-12-31	0.027067	3.954405	105.571953
1852-12-31	0.028324	3.955504	105.392555
1853-12-31	0.028730	3.955512	106.455231
1854-12-31	0.027566	3.956528	109.545181

Visualizzazione

# Costruzione di un dataset standardizzato per il plot
from sklearn.preprocessing import StandardScaler
ts_zosStd = pd.DataFrame(StandardScaler().fit_transform(ts_zos), columns=ts_zos.columns)
ts_zosStd = ts_zosStd.set_index(ts_zos.index)

plt.rcParams.update({'figure.figsize': (9, 5)})
ts_zosStd.plot()
plt.title("Confronto tra le serie storiche standardizzate delle tre variabili")
plt.show()

Standardizzando le tre variabili rendiamo possibile un confronto tramite una rappresentazione su grafico.
Le variabili "zos" (livello del mare rispetto al geoide) e "bigthetaoga" (temperatura conservativa del mare) sono visibilmente molto simili. Entrambe infatti seguono un trend positivo molto simile nella prima metà dello span temporale. Una differenza però si dimostra dagli anni 1960: mentre il trend della temperatura conservativa rimane lineare, quello della variabile di riferimento accellera, diventando quasi esponenziale.
La variabile relativa alla massa dei ghiacciai invece si dimostra la più diversa visivamente. Ciò è dovuto a, come accennato prima, la sua presunta correlazione negativa rispetto alle altre.

Calcolo della correlazione tra le variabili

ts_zos.corr()

	zos	bigthetaoga	simass
zos	1.000000	0.937414	-0.495248
bigthetaoga	0.937414	1.000000	-0.638905
simass	-0.495248	-0.638905	1.000000

La matrice di correlazione porta buoni risultati: la nostra variabile di riferimento "zos" è correlata al 93.74% con la variabile "bigthetaoha" e al -49.52% con la variabile "simass". Possiamo visualizzare le informazioni contenute nella matrice tramite una heat-map:

# Build the corplot
plt.imshow(ts_zos.corr(), cmap = 'RdBu_r', interpolation='nearest', vmin=-1)
# Add legend
plt.colorbar()
# Create list with tick_mark positions
tick_marks = [i for i in range(len(ts_zos.columns))]
# Add the tickmarks at the designated position, using the labels of the dataframe
plt.xticks(tick_marks, ts_zos.columns, rotation='vertical')
plt.yticks(tick_marks, ts_zos.columns)
plt.title("Heat-Map della matrice delle correlazioni tra i tre modelli")
plt.show()

Granger Causality test

A questo punto, valutiamo la capacità previsiva delle serie storiche da noi costruite, per fare ciò utilizziamo il test di causalità di Granger. La causalità di Granger è un concetto espresso nel 1969 da Clive Granger e ampliato successivamente da Christopher Sims mirante a determinare una causalità tra le variabili espresse in un modello. Una serie storica $\{x_t\}_t$ causa (nel senso di Granger) una seconda serie storica $\{y_t\}_t$ se, condizionando rispetto ai valori passati di $x_t$, l’errore quadratico medio di previsione della $y+t$ risulta ridotto rispetto al caso in cui l’informazione dei valori passati di $𝑥𝑡$ sia ignorata. Formalmente: $$ E[y_t − E(y_t| \cdot) | y_{t-1}, y_{t−2}, … ; x_{t-1}, x_{t-2}, … ]^2 \le E[y_t − E(y_t| \cdot) | y_{t−1}, y_{t−2}, … ]^2 $$

from statsmodels.tsa.stattools import grangercausalitytests

test = grangercausalitytests(ts_zos[['zos','bigthetaoga']], 30, verbose=False)
for i in test:
    if test[i][0]['ssr_ftest'][1] < 0.05:
        print('found significant p-value at lag ' + str(i) + ': ' + str(round(test[i][0]['ssr_ftest'][1] , 3)))

test = grangercausalitytests(ts_zos[['zos','simass']], 30, verbose=False)
for i in test:
    if test[i][0]['ssr_ftest'][1] < 0.05:
        print('found significant p-value at lag ' + str(i) + ': ' + str(round(test[i][0]['ssr_ftest'][1] , 3)))

Il test di granger, però, non trova alcun valore significativo a nessun livello di ritardo per nessuna delle due variabili. Per questo motivo, nonostante le serie siano molto correlate, non possiamo dire che né la temperatura conservativa né l'area dei ghiacciai siano utili a prevedere l'innalzamento dei mari.

Aggiunta di un ultimo modello

Per determinare se gli scarsi risultati avuti dal Granger Test siano legati o meno al modello scelto abbiamo deciso di ripetere le analisi sostituendo il modello del "livello del mare rispetto al geoide" con la variabile "altezza del mare globale osservata".
Per fare ciò, utilizziamo il modello presente al seguente link sviluppato durante il seguente progetto: "Frederikse, T., Landerer, F., Caron, L. et al. The causes of sea-level rise since 1900. Nature 584, 393–397 (2020)" link

# Apro ed esploro il modello 4
modfile4 = 'C:/Users/david/Desktop/GeoData/Project/global_timeseries_measures.nc'
d4 = xr.open_dataset(modfile4)
d4

<xarray.Dataset>
Dimensions:                                             (time: 119)
Coordinates:
  * time                                                (time) datetime64[ns] ...
Data variables: (12/21)
    global_average_sea_level_change                     (time) float32 -173.3...
    global_average_sea_level_change_upper               (time) float32 -141.0...
    global_average_sea_level_change_lower               (time) float32 -203.6...
    glac_mean                                           (time) float32 -79.55...
    glac_upper                                          (time) float32 -62.5 ...
    glac_lower                                          (time) float32 -95.5 ...
    ...                                                  ...
    global_average_thermosteric_sea_level_change        (time) float32 -48.28...
    global_average_thermosteric_sea_level_change_upper  (time) float32 -37.55...
    global_average_thermosteric_sea_level_change_lower  (time) float32 -58.67...
    sum_of_contrib_processes_mean                       (time) float32 -173.1...
    sum_of_contrib_processes_upper                      (time) float32 -146.8...
    sum_of_contrib_processes_lower                      (time) float32 -199.2...
Attributes: (12/42)
    title:                     Global sea-level changes and contributing proc...
    summary:                   This file contains reconstructed global-mean s...
    id:                        10.5067/GMSLT-FJPL1
    naming_authority:          gov.nasa.jpl
    source:                    Frederikse et al. The causes of sea-level rise...
    project:                   NASA sea-level change science team (N-SLCT)
    ...                        ...
    time_coverage_start:       1900-01-01
    time_coverage_end:         2018-12-31
    time_coverage_duration:    P119Y
    time_coverage_resolution:  P1Y
    date_created:              2020-07-28
    date_modified:             2020-09-14

plt.rcParams.update({'figure.figsize': (9, 5)})
d4.global_average_sea_level_change.plot()
plt.show()

Il trend della variabile non si discosta troppo da quello della variabile usata in precedenza. Ci accorgiamo però che presenta un andamento meno lineare, soprattutto nell'ultima parte in cui sembrerebbe accellerare, proprio come la variabile relativa alla temperatura conservativa del mare.

# Rimozione warning
pd.options.mode.chained_assignment = None

# Costruzione di un dataset che contenga tutte le variabili
ts_zos2 = ts_zos.drop('zos', axis=1)
ts_zos2 = ts_zos2[50:]
ts_zos2['zos'] = d4.global_average_sea_level_change.values.tolist()[:-4]
# Cambio dell'ordine delle colonne
ts_zos2 = ts_zos2[['zos', 'bigthetaoga', 'simass']]
ts_zos2.head()

	zos	bigthetaoga	simass
year
1900-12-31	-173.261398	3.975367	108.662781
1901-12-31	-168.801697	3.975936	109.415802
1902-12-31	-180.389999	3.975817	117.218178
1903-12-31	-170.267807	3.974838	117.983482
1904-12-31	-168.720093	3.974608	113.168900

Visualizzazione

# Costruzione dataset standardizzato per visualizzazione
ts_zosStd2 = pd.DataFrame(StandardScaler().fit_transform(ts_zos2), columns=ts_zos2.columns)
ts_zosStd2 = ts_zosStd2.set_index(ts_zos2.index)

plt.rcParams.update({'figure.figsize': (9, 5)})
ts_zosStd2.plot()
plt.title("Confronto tra le serie storiche standardizzate delle tre variabili")
plt.show()

Calcolo delle Correlazioni tra le variabili

ts_zos2.corr()

	zos	bigthetaoga	simass
zos	1.000000	0.918667	-0.777208
bigthetaoga	0.918667	1.000000	-0.849646
simass	-0.777208	-0.849646	1.000000

# Build the corplot
plt.imshow(ts_zos2.corr(), cmap = 'RdBu_r', interpolation='nearest', vmin=-1)
# Add legend
plt.colorbar()
# Create list with tick_mark positions
tick_marks = [i for i in range(len(ts_zos2.columns))]
# Add the tickmarks at the designated position, using the labels of the dataframe
plt.xticks(tick_marks, ts_zos2.columns, rotation='vertical')
plt.yticks(tick_marks, ts_zos2.columns)
plt.title("Heat-Map della matrice delle correlazioni tra i tre modelli")
plt.show()

Granger Causality Test

test = grangercausalitytests(ts_zos2[['zos','bigthetaoga']], 30, verbose=False)
for i in test:
    if test[i][0]['ssr_ftest'][1] < 0.05:
        print('found significant p-value at lag ' + str(i) + ': ' + str(round(test[i][0]['ssr_ftest'][1] , 3)))

found significant p-value at lag 1: 0.01
found significant p-value at lag 12: 0.048
found significant p-value at lag 14: 0.033

test = grangercausalitytests(ts_zos2[['zos','simass']], 30, verbose=False)
for i in test:
    if test[i][0]['ssr_ftest'][1] < 0.05:
        print('found significant p-value at lag ' + str(i) + ': ' + str(round(test[i][0]['ssr_ftest'][1] , 3)))

found significant p-value at lag 20: 0.046
found significant p-value at lag 21: 0.024
found significant p-value at lag 22: 0.025
found significant p-value at lag 23: 0.039

I risultati del Granger Causality Test ottenuti utilizzando la nuova variabile dipendente sono decisamente più interessanti.
Il test trova infatti p-value significativi a ritardi 1, 12 e 14 tra le variabili "zos" e "bigthetaoga". Ciò significa che la temperatura conservativa dei mari è utile a predirre l'innalzamento dei mari a ritardi di 1, 12 e 14 anni.
Per quanto riguarda la variabile "simass", invece, i risultati più interessanti si mostrano a ritardi più elevati. C'è infatti evidenza di una capacità della variabile relativa alla massa di ghiaccio nel mare di predirre l'innalzamento del livello dei mari a ritardi variabili tra 20 e 23 anni.